rm(list=ls())

library('reshape2')
library('dplyr')
library('matrixStats')
library('ggplot2')



CFpe <- read.csv("matlab/baseline_model/selection_prpeace_replication.csv")
CFpe[,3:6]  <- 100*CFpe[,3:6]

cf_labs <- c("Max peace", "Max conflict", "Best rebel")
cv <- 1.96

Mpe <- as.matrix(select(CFpe, -one_of(c("ccode","countrycode","data"))))



ggdat <- data.frame(means = colMeans(Mpe),
                     sds = colSds(Mpe),
                     cf = cf_labs,
                     var = rep("Probability of Peace", each=3))
ggdat2 <- data.frame(var="Probability of Peace",
                    Z=c(mean(CFpe$data)),
                    W = c(sd(CFpe$data)),
                    X = 0.5,
                    Y = 0.75)

ggdat_use <- subset(ggdat, 
                    cf %in% c("Max peace", "Max conflict") &
                    var == "Probability of Peace")

plot_peace <- ggplot(ggdat_use) + 
  geom_pointrange(aes(x=factor(cf, levels=cf_labs), y=means, ymax=means+cv*sds, ymin=means-cv*sds), size=0.75) + 
  theme_bw(18) + ylab("Probability of Peace") + xlab("Selection Mechanism")  + 
  geom_hline(data=ggdat2[1,], aes(yintercept = Z), linetype="dashed", color="black", size=1.2) + 
  geom_rect(data = ggdat2[1,], aes(ymin = Z-cv*W, ymax = Z+cv*W),
            xmin=1-1,
            xmax=3+1,
            fill="black", alpha=0.25) +
  theme(axis.title.x=element_text(margin = margin(t = 15, r = 0, b = 0, l = 0))) 

ggsave("fig3.pdf", plot_peace, height=4, width=6.5)